program powerjb
use moments3jb
use matrixfuns
use datahadler
use linear_operators
use imhoff
implicit none
integer, parameter:: n=3,p=6,num=100,T=5000,pth=9,s=10000
double precision,parameter:: si=.999d0
double precision wmat(p,p),eta(num),power(num,3)& !joint, skewness,kurtosis
&,pr,crit,critk,critsk,dchiin&
&,b(n),iwmat(p,p),matm(p,1),matv(p,p)
double precision sigma(n,n),MD(p+pth),VD(p+pth,p+pth)&
&,dmut(pth,n),dvsig(pth,n*n),dmat(n,n),matA(pth,pth),matB(p,pth),matD(p,pth+p)&
&,powers(num),rej
integer i,j,k,ifault,irank



forall(i=1:pth,j=1:n)
dmut(i,j)=0.0d0
end forall
forall(i=1:pth,j=1:n*n)
dvsig(i,j)=0.0d0
end forall
forall(i=1:n)
	dmut(i,i)=1.0d0
end forall
forall(i=1:n,j=1:n)
	dmat(i,j)=0.0d0
end forall
k=0
do j=1,n,1
	do i=1,n,1
		k=k+1
		dmat(j,i)=1.0d0
		dmat(i,j)=1.0d0
		dvsig(n+1:pth,k)=vech(dmat)
		dmat(j,i)=0.0d0
		dmat(i,j)=0.0d0
	end do
end do


eta=linspace(.0001d0,.02d0,num)
b=vassign(n,0.75d0)

sigma=eye(3)
!sigma(1,:)=(/1.0d0,0.1d0,0.2d0/)
!sigma(2,:)=(/0.1d0,1.0d0,0.3d0/)
!sigma(3,:)=(/0.2d0,0.3d0,1.d0/)


crit=dchiin(.95d0,dble(p))
critk=dchiin(.95d0,dble(n))
critsk=dchiin(.95d0,dble(n))

call var0u(wmat,dmut,dvsig,pth,sigma,b)
iwmat=.i.wmat

do i=1,num,1
	call normjbmvu(MD,VD,matA,matB,b,eta(i),si,sigma,dmut,dvsig,pth)
	matD(:,pth+1:pth+p)=eye(p)
	matD(:,1:pth)=matmul(matB,.i.matA)

	matm(:,1)=dsqrt(dble(T))*matmul(matD,MD)
	matv=(matD.x.VD).xt.matD
!Joint
	call imhofprob(pr,ifault,matm,matv,iwmat,crit)
	power(i,1)=1.0d0-pr
if (ifault .ne. 0) then
	print*, "joint",ifault
end if
!Skewness
	call imhofprob(pr,ifault,matm(n+1:p,:),matv(n+1:p,n+1:p),.i.wmat(n+1:p,n+1:p),critsk)
	power(i,2)=1.0d0-pr
if (ifault .ne. 0) then
	print*, "skewness",ifault
end if
!Kurtosis
	call imhofprob(pr,ifault,matm(1:n,:),matv(1:n,1:n),.i.wmat(1:n,1:n),critk)
	power(i,3)=1.0d0-pr
	
	print*, i,real(power(i,:))
if (ifault .ne. 0) then
	print*, "kurtosis",ifault
end if
end do
call exportvec(eta,"ejb.txt",1)
call exportvec(vec2(power),"powrjb3a.txt",1)

end program powerjb